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ABSTRACT 


Nine  explosions  have  been  used  to  study  the  problem  of  specifying  the  num¬ 
ber  and  characteristics  of  the  arrival  phases  of  a  nuclear  explosion.  Tenta¬ 
tive  identification  of  the  £ui^litudes  and  time  delays  of  the  multiple  arrivals 
are  given. 

The  method  used  is  to  scan  the  likelihood  (L)  for  a  single  reflection  model 
over  eunplitude  (a)  and  time  delay  (t ) .  The  likelihood  is  constructed  in  the 
frequency  domain  using  von  Seggern  and  Blandford's  (1972)  source  signal  spectrum 
and  a  noise  spectrum  estimated  from  a  time  window  preceding  the  signal. 

Ihe  results  denonstrated  that  pP  was  present  along  with  additional  multiple 
arrivals  indicated  by  other  peaks  in  the  likelihood  function,  L.  (Analysis 
extended  to  delays  of  1.5  seconds,  too  short  to  detect  spall  except  for  PILE- 
DRIVER,  for  which  case  it  could  not  be  detected.)  In  the  case  of  MAST  and 
PILEDRIVER,  the  likelihood  estimate  for  x  was  found  to  be  clearly  superior  to 
T  «  l/^n*  where  f^  is  the  first  null  frequency.  The  result  for  PILEDRIVER, 

T  “  .15,  is  the  first  accurate  estimate  of  x  for  this  event,  £md  the  shortest 
estimated  delay  time  in  the  literature.  The  delay  time  for  CANNIKIN,  .55  sec, 
is  significantly  shorter  than  twice  the  observed  uphole  time:  1.0  sec.  This 
result  is  unexplained,  but  it  may  be  related  to  a  change  of  reflection  coef¬ 
ficient  with  frequency  which  is  not  included  in  the  model.  Con$>arison  of 
theoretical  and  observed  WHSSN  long-period  explosion  P  waves  provides  evidence 
for  such  a  variable  reflection  coefficient. 

Fair  to  good  agreement  was  obtained  between  predicted  euid  observed  explo¬ 
sion  periods  measured  on  WWSSN  long-period  film  and  on  LRSM  short-period  data 
which  was  deconvolved  to  appear  as  if  recorded  on  the  WWSSN  long-period  system. 

A  promising  technique  to  determine  yield  and  pP  delay  time  given  the  WWSSN 
long-period  eunplitude  Md  period  is  presented. 

All  these  approaches  require  that  both  t*  =  travel-time/Q,  and  the  reduced 
displacement  potential  as  a  function  of  depth  and  yield  be  known. 
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INTRODUCTION 


0«tennlning  the  P  -  pP  delay  poses  an  interesting  problem  during  analy¬ 
sis  of  data  generated  by  nuclear  explosions.  This  delay  can  be  used  to  pro¬ 
vide  an  estimate  for  the  depth  of  the  event. 


The  model  of  interest  is  that  the  P-wave  signal  will  be  seen  on  the 


seismometer  as 


St  ^  as^ 


+  n 
T  t 


(1.1) 


where  t  is  the  delay  time  for  pP  and  £  is  the  relative  amplitude  for  the 
reflected  wave.  The  problems  of  estimating  amplitude  and  the  delay  time 
usually  depends  upon  noticing  that  for  stationary  and  independently  distri¬ 
buted  signal  and  noise  processes  with  spectra  P  (f)  and  P  (f) ,  respectively, 

s  n 

the  spectrum  of  observed  data  in  (1.1)  is  given  as 

2. 


where 


P  (f)  =  |A(f) I  P  (f)  +  P  (f) 
y  s  n 


A(f)  =  (1  +  ae"^’'  '^) 


Here,  we  would  have 


|A(f)  1^  =  1  +  a^  +  2a  cos  (2TrfT) 


(1.2) 


(1.3) 


(1.4) 


For  a  =  -1  in  Equation  (1.4),  as  might  be  expected  for  pP,P^(f)  will  be 


zero  for 


f  =  -,  m  =  0,  1, , 

T 


Cohen  (1970) ,  using  detection  of  nulls,  estimated  the  time  delay  as  the  re¬ 
ciprocal  of  the  frequency  for  the  first  spectral  null.  Figure  1  shows  the 
spectrum  without  nulls.  Figure  2  shows  the  spectrum  of  the  artificial  time 
series 

and  we  observe  a  first  null  at  1.56  Hz  so  that  .64  seconds  is  the  estimated 
delay  time.  The  distortion  due  to  the  shape  of  the  theoretical  power  spectrum 
of  the  signal  P^(f),  produced  the  discrepancy  between  0.64  and  0.7  in  this 


case. 


The  normalized  or  "whitened"  spectrum  in  which  Pg(f)  is  divided  out 


showed  a  local  minimum  at  .70  seconds.  Cohen  (1970)  used  this  whitening 
method  to  determine  pP  and  spall-created  secondary  phase  delays. 


Cohen,  T.  J.  (1970),  So\jrce-depth  determination  using  spectral  pseudo- 
autocorrelation  and  cepstral  analysis,  Geophys.  j.  R.  Astr.  Soc.. 
20,  223-231.  - 


Another  method  of  analysis  (Bogert  et  al. ,  1962)  is  based  upon  the  ob¬ 
servation  that  for  the  noiseless  case,  the  logarithm  of  the  spectrum  (1,2) 
is  approximately,  if  a  is  small, 

2 

logP  (f)  =  logP  (f)  +  logd  +  a  +  2acos  2  irfr) 
y  s 

2 

=  logP  (f)  +  a  +  2acos2  wfx  (1.5) 

s 

implying  that  the  log  spectrum  is  periodic.  Thus,  the  spectrum  of  logP  (f) 
should  peak  at  t  =  1/f.  Ideally  one  would  detrend  by  subtracting  the  log  of 
the  source  spectrum  as  well  as  the  overall  mean.  Cohen  (1970)  and  Flinn  et 
al.  (1973)  observed  that  for  larger  echo  amplitudes  (|a|  >  .5)  the  spectrum 
of  the  original  spectrum  produces  a  more  pronounced  peak  than  the  spectrum 
of  the  log  spectrum. 


Bakun  and  Johnson  (1973)  fitted  spectra  (1.2)  using  least  squares  tech¬ 
niques  to  the  CANNIKIN  and  MILROW  events,  identifying  pP  and  spall  phases 
along  with  the  respective  aunplitudes  and  delay  times.  They  also  did  the 
analysis  in  the  time  domain  using  homomorphic  deconvolution  to  produce  pulse 
wavelets  corresponding  to  the  different  phases. 


A  natural  extension  of  least  squares  fitting  is  to  apply  maximum  like¬ 
lihood  in  the  frequency  domain  to  develop  estimators  for  the  time  delays  and 
the  amplitudes.  Hannan  and  Thomson  (1974)  wrote  the  likelihood  function 
assuming  that  the  delay  times  and  amplitudes  are  frequency  dependent. 


Bakun,  W.  H.  and  L.  R.  Johnson  (1973).  The  deconvolution  of  teleseismic  P 
waves  from  explosions  MILROW  and  CANNIKIN,  leophys.  J.  R.  Astr.  Soc., 

34,  321-342. 

Bogert,  B.  P.,  Healy,  M.  J.  R.  and  Tukey,  J.  W.  (1962).  The  frequency 

analysis  of  time  series  for  echoes:  cepstrum,  pseudcf-autocovariance , 
cross  sepstrum  and  phase  cracking,  in  Proc.  on  Time  Series  Analysis, 
ed.  M.  Rosenblatt,  John  Wiley  and  Sons,  New  York. 

Flinn,  E.  A,,  Cohen,  T.  J.  and  McCowan,  D.  W.  (1973).  Detection  and  analysis 
of  multiple  seismic  events.  Bull.  Seis.  Soc.  Am., 63 (6),  1921-1936. 


Hannan,  E.  J.  and  Thomson,  P.  J.  (1974) .  Estimating  echo  times.  Technometrics 
16(1) ,  p.  77-84. 


However,  they  observed  that  frequency  dependent  estimators  for  these  parame¬ 
ters  tended  to  be  about  the  same  for  each  frequency.  Estimators  for  the  siq- 
nal  eind  noise  spectra  are  made  over  narrow  bands  assuming  that  they  are  con- 
stant  over  the  narrow  bands  and  that  |A(f)|  varies  significantly  over  the 
same  band.  I!owever,  these  assumptions  are  not  reasonable  for  the  seismic 
data  considered  here  because  in  this  study  the  authors  assume  in  most  cases 
that  the  source  signal  spectrum  follows  the  model  that  von  Seggern  and 
Blandford  (1972)  described  and  that  an  independent  estimator  for  noise  spec¬ 
trum  can  be  derived  from  a  time  window  preceding  the  signal.  The  approximate 
likelihood,  then,  can  be  written  in  terms  of  the  delays  and  amplitudes  over 
the  whole  frequency  band  of  interest. 

The  basic  data  analyzed  consists  of  nine  selected  explosions  as 
recorded  at  several  LRSM  stations  (Table  I).  The  spectrum  analyzed  is 
determined  as  follows;  Good  signal-to-noise  traces  are  selected.  A  cosine 
window  is  applied  to  the  two  seconds  of  noise  preceding  first  motion,  and 
another  cosine  taper  is  applied  to  the  last  five  seconds  of  the  records. 

Each  of  the  traces  is  12.8  seconds  long.  A  12.8  second  noise  sample  in 
in  front  is  identically  tapered.  The  signal  traces,  preceded  by  their 
associated  noise,  are  normalized  by  the  maximum  signal  amplitude.  To  obtain 
the  final  "network"  spectrum,  the  median  of  the  signal  spectra  and  of  the 
noise  spectra  is  calculated  at  each  frequency. 


von  Seggern,  D.  and  Blandford,  R.  (1972).  Source  time  functions  and 
for  underground  nuclear  explosions,  Geophys.  J.  R.  A.  Soc.,  31, 


sp>ectra 

83-97. 
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Faraaeters  of  explosions  considered  in  this  paper,  together  with 
LRSM  seisao^ras  nusibers  and  start  points  for  the  data  ised 


LIKELIHOOD  FUNCTION 


The  model  for  the  time  series  can  be  written  as  (1.1),  leading  to  a 
model  of  the  form  (1.2),  (1.3)  for  the  spectrvim.  Suppose  that  an  arbitrary 
signal  gain  factor,  0  is  incorporated  into  the  original  equation  yielding 


and 

\  = 

e(St  + 

(2.1) 

Py(f) 

=  9^  |A(f) P  (f)  +  P  (f) 
s  n 

(2.2) 

where 

lA(f) 

1  2  1  -2  TTif  T  1  2 

1  =  1  1  +  ae  1 

2 

=  1  +  a  +  2acos  (2  ir  f  t  ) 

(2.3) 

The  signal  spectrum  P  (f)  is  derived  from  the  von  Seggern  and  Blandford 

s 

(1972)  model  modified  for  the  effects  of  the  earth's  ^ibsorption,  that  is 


.f...  - 


(1  +  (1+2B)^  ( 


.)2)l/2 


e  ~  ’^^^*I(f) 


(2.4) 


(1  +  (^-^)^) 


where  1(f)  is  the  instrument  response,  t*  =  T/Q,  where  T  is  the  travel-time, 
and  Q  the  quality  factor.  The  parameter  C  is  related  to  yield  Y  in  kilotons 
through  the  equation 

C  =  k^(y/5)"^^^  (2.5) 

The  observed  spectrum  is,  then,  a  function  of  the  parameters  in  the  signal 

2 

spectrum  B,  C,  t*,  as  well  as  the  gain  0  and  amplitude-time  delay  pair  a,  t. 

While  B  is  a  non-dimensional  parameter,  C,  k  have  the  units  of  time. 

—  o 

Now,  let  Y^(f)  be  the  FFT 


y  (f)  =  Ey 

t=0  '  (2.6) 

so  that  an  approximation  to  the  joint  log  likelihood  over  L  adjacent 
frequencies  F  which  are  integer  multiples  of  -  can  be  written  as 

von  Seggern,  D.  and  Blandford,  R.  (1972) .  Source  time  functions  and 

spectra  for  underground  nuclear  explosions,  Geophys.  J.  R.  Astr.  Soc. 
31,  83-97. 
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L(e  ,  a,  t,  B,  c,  t*) 


=  -Llogti  -  ElogP  (f) 
f 


E 

f  PyCf) 


Assume  first  that  parameters  B,  C,  and  t*  can  be  determined  reasonably  well 

from  site  characteristics.  An  initial  approximation  to  the  maximum  likeli- 

2 

hood  estimator  for  0  is  derived  by  setting  the  derivative  of  (2.7)  equal  to 
0  and  assuming  that  '  0-  This  procedure  yields 

52.  1e  I' 

'  |A(f)  |\(f) 

(2.8) 

where  the  solution  depends  upon  the  assvmed  values  of  a  and  t.  The 

^2 

likelihood  (2.7)  is  then  scanned,  with  9  changing,  over  the  values  of 
a  and  T  restricted  to  a  reasonable  set,  such  as  -1  <  a  <  1,  .05  sec 
^  T  ^2  sec.  The  range  of  f  in  (2.7)  and  (2.8)  is  0.3  to  2.0  Hz  in 
order  to  insure  that  only  data  with  high  values  of  S/N  are  considered. 

Note  that  |v^(f) 1  ^  is  just  the  usual  squared  FFT  or  the  sample  periodogram. 

In  order  to  illustrate  use  of  the  likelihood  function,  artificial 
data  was  generated,  conforming  to  the  models 


Yt  -  2S^  + 


(2,10) 


V  ■  -  -“t  -  .7  *  -  .25  -  -  .95]  >  *  "t 

While  the  first  equation  has  only  a  single  arrival,  (2.10)  has  a  pP 
arrival  at  t  =  0.7  seconds.  Equation  (2.11)  has  a  pP  arrival  at  0.7 
seconds  and  then  a  multipath  of  the  complete  signal  .25  seconds  later. 
In  all  cases  we  used  B(=2.0,  C=7.6  (this  corresponds  to  Y=10  kilotons  in 


4 


granite)  and  t*  =  .55  as  known  parameters  in  the  signal  model  (2.4).  The 
theoretical  process  white  and  negligible. 

Figure  1  shows  the  likelihood  function  under  (2.9);  note  that  it  attains 
a  maximum  corresponding  to  a=0  in  model  (2.1)  for  any  time  delay,  which  leads 
to  the  single  arrival  explanation. 

When  the  single  reflection  model  (2.10)  is  introduced,  the  likelihood 
in  Figure  2  shows  a  reasonably  well  defined  peak  at  a  =  -.5  and  t  =  .7  seconds. 
One  observed  characteristic  is  that  the  time  delay  tends  to  be  determined 
with  higher  resolution  than  the  amplitude.  Figure  2  shows  that  the  observed 
and  fitted  spectrum  agree  very  well. 

The  likelihood  function  even  though  it  was  designed  for  a  single  re¬ 
flection  model,  can  easily  handle  the  more  complicated  double  arrival  defined 
by  Equation  (2.11).  Figure  3  demonstrates  that  the  theoretical  spectrum  of 
the  data  (P^Cf))  shows  no  rippling  for  this  signal  spectrum.  Thus,  neither 
spectral  analysis  nor  looking  for  nulls  is  expected  to  do  well  in  this  situa¬ 
tion.  Nevertheless,  the  two  peaks  in  the  likelihood  correspond  to  the 
correct  amplitudes  and  delay  times.  A  narrowing  of  the  main  peak  is  the  main 
impact  on  the  signal  spectrum. 

The  simulated  examples  seem  to  confirm  that  if  the  estimated  source 
spectrum  is  correctly  specified,  then  the  multiple  reflections  can  be  sorted 
out  even  though  they  may  not  introduce  any  discernible  periodicity  in  the 
observed  spectrian.  While  the  problem  of  independently  estimating  the  B,  C, 
and  t*  parameters,  using  the  likelihood  function,  has  not  been  considered  in 
this  initial  study,  it  will  be  investigated  in  the  following  section  for  a 
number  of  nuclear  explosions. 
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Figure  1.  Analysis  of  artificial  data  generated  according  to 
the  model  y.  =  2S  +  n  . 
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Figure  2.  Analysis  of  artificial  data  generated  according  to 
the  model  y  =  2(S  -  .5S  )  +  n  . 
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Figure  3.  Analysis  of  artificial  data  generated  according  to 
the  model  =  2  (S^  -  .5S^  ^  +  .5(S.  .j.  -  .5S 


ESTIMATION  OP  pP  DELAYS  AND  DEPTHS  FOR  NINE  EXPLOSIONS  USING  LIKELIHOOD 


Table  I  shows  the  nine  explosions  analyzed  using  the  likelihood  function 
along  with  information  on  depth  and  velocity  in  the  overburden.  The  observed 
spectra  are  the  median  spectra  from  various  combinations  of  LRSM  stations  as 
discussed  earlier.  Figures  4  through  12  display  the  likelihood  functions,  the 
observed  spectra,  and  the  fitted  spectra  corresponding  to  combinations  of  ar¬ 
rivals  appropriate  to  sets  of  the  observed  maxima. 

Table  II  summarizes  estimated  arrivals  for  each  event  with  the  arrivals 
selected  as  pP  denoted  by  *.  In  most  cases,  pP  is  identified  by  looking  for 
an  arrival  maximizing  the  likelihood  and  appearing  with  a  negative  amplitude. 

For  example,  the  peaks  in  the  likelihood  function  indicate  four  arrivals 
corresponding  to  MILROW  (Figure  4)  with  the  third  arrival  maximizing  the  like¬ 
lihood  function  at  an  amplitude  of  -.45  and  time  delay  of  .85  seconds.  This 
delay  is  consistent  with  the  first  null  observed  at  .80  seconds.  An  early 
poorly  understood  reflection  (multipath?),  which  appeared  at  .25  seconds  with 
a  positive  amplitude,  showed  frequently  with  these  events.  Theoretical  spectra 
in  Figure  4  are  synthesized  assuming  two  or  four  arrivals  with  the  delay  times 
and  amplitudes  indicated.  The  two-arrival  model  appears  to  yield  a  better 
fit  for  this  particular  case.  In  principal,  each  positive  amplitude  arrival 
(multipath)  could  have  a  corresponding  negative  amplitude  for  pP.  However, 
these  eunplitudes  would  be  small  and  they  have  been  neglected  in  the  recon¬ 
structions  . 

In  each  case  except  BUKARA  which  was  detonated  in  clay  (Marshall,  1972), 
the  scaled  granite  source  model  of  von  Seggern  and  Blandford  was  used.  An 
atten^t  to  match  the  overall  spectral  shape  of  BUKARA  by  visual  matching  of 
observed  and  theoretical  si^ectra  led  to  the  source  spectral  values  B=0,  C=4.0. 
Matching  the  high  frequency  (2-5  Hz)  slopes  of  the  median  spectra  with  good 
signal-to-noise  ratios  led  to  t*  values  of  .55  for  NTS  explosions,  and  .27 
for  Amchitka  explosions.  These  should  be  regarded  as  "median"  t*  values 
suitable  for  the  "median"  spectrum. 


Marshall,  P.  D.,  (1972).  Some  seismic  results  from  a  world  wide  sample  of 
large  underground  explosions,  AWRE  Report  0  49/72,  Atomic  Weapons 
Research  Establishment,  Aldermaston,  Berkshire,  United  Kingdom. 
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Figure  4.  Analysis  of  MILROW  event  showing  pP  at  0.85  sec 
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Figure  7.  Analysis  of  LONGSHOT  event  showing  pP  at  0.55  sec 
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gure  9.  Analysis  of  EDAM  event  shoving  pP  at  0.55  sec 
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10.  Analysis  of  MAST  event  showing 
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Analysis 


TABLE  11 


Parameters  of  maximum-likelihood  estimated  arrival  phases  for  a  suite 
of  explosions,  together  with  the  B,  C,  and  t*  values  assumed 


a 

T 

L 

B 

C 

1. 

MILROW 

1 

.40 

.30 

1.377 

2.00 

2.87 

Estimated  from 

2 

-.30 

.55 

1.369 

first  null  t  = 

*3 

-.45 

.85 

1.378 

(.80) 

4 

.35 

1.40 

1.374 

2. 

CANNIigN 

1 

.35 

.30 

1.428 

2.00 

1.68 

Estimated  from 

*2 

-.35 

.55 

1.429 

first  null  r  = 

3 

.15 

1.70 

1.415 

(1.16) 

4 

-.25 

2.05 

1.422 

5 

-.25 

2.25 

1.422 

3. 

GREELEY 

1 

.60 

.25 

1.771 

2.00 

2.87 

No  visible 

2 

-.45 

.50 

1.701 

null 

*3 

-.55 

.70 

1.712 

4. 

LONGSHOT 

1 

.50 

.25 

1.440 

2.00 

6.67 

Estimated  from 
first  null  t  = 
(.51) 

*2 

-.60 

.55 

1.450 

5. 

BUKARA 

1 

.55 

.35 

0.00 

4.00 

Estimated  from 
first  null  t  = 
(1.07) 

*2 

-.75 

1.20 

1.314 

6. 

EDAM 

*1 

-.40 

.55 

1.487 

2.00 

(k  =16.8) 

Estimated  from 
first  null  t  = 
(.67) 

2 

.45 

1.20 

1.486 

o 

7. 

MAST 

1 

.20 

.25 

1.551 

2.00 

(k  =16.8) 

Estimated  from 

*2 

-.35 

.50 

1.564 

o 

first  null  t  = 

3 

-.30 

.95 

1.554 

(.80) 

4 

.30 

1.30 

1.562 

5 

-.10 

1.65 

1.549 

8. 

MIZZEN 

1 

.45 

.30 

1.463 

2.00 

(k  =16.8) 

Estimated  from 
first  null t  = 
(.91) 

*2 

-.60 

.95 

1.453 

o 

9. 

PILEDRIVER 

*1 

-.5 

.15 

1.437 

2.00 

7.60 

Estimated  from 
first  nullT  = 
(.64)7 

2 

.3 

.'40 

1.434 

* 

Phase 

selected  for 

pP 

-  27  - 

-  1 

1 

.27 


.27 


.55 


.27 


.55 


.55 


.55 


A  computation  for  the  apparent  depth  can  be  made  using  d  =  v  t/2  or  for 
the  predicted  delay  by  t  =  2(d/v)  where  v  is  the  shot-po'nt  to  surface  average 
velocity  of  the  overburden  (Table  I)  and  t  is  the  estimated  pP  time  delay. 
Figure  13  shows  the  depths  computed  from  the  likelihood  approach  compared  v;itl. 
published  depths  for  the  nine  events. 

Below,  the  events  are  discussed  separately; 

MILPOM  (Figure  4) 

The  MILROW  results  are  in  good  agreement  with  Bakun  et  al's  (1973) 
earlier  results  and  with  the  pP  estimates  from  Springer's  (1974)  shot  point 
recordings.  Depths  estimated  from  the  first  null  and  from  maximum  likelihood 
are  comparable. 

CANNIKIN  (Figure  5) 

For  CANNIKIN  achieving  good  agreement  is  possible  with  the  observed 
spectrum  with  a  suite  of  five  arrivals,  including  one  with  negative  amplitude 
and  with  maximiun  likelihood  at  0.55  seconds  delay.  Springer  (1974)  obtained 
1.0  seconds  for  pP,  and  Bakun  and  Johnson  (1973)  obtained  1.5  seconds.  Con¬ 
sideration  of  the  first  null  leads  to  an  estimated  delay  of  1.16  seconds. 
Clearly,  our  result  disagrees  with  previous  work.  Use  of  delays  and  amplitudes 
proposed  by  Bakun  and  Johnson,  which  include  a  multipath,  results  in  a 
theoretical  spectrum  that  is  not  as  good  a  match  to  our  observed  spectraum  as 
the  theoretical  spectrum  determined  in  this  report.  Adjustments  of  t*  have 
been  attempted ,  but  they  did  not  yield  a  maximum  likelihood  estimate  for  a 
negative  delay  near  1.0  second.  This  result  suggests  a  failure  of  the  maximxim- 
likelihood  technique.  One  explanation  for  this  outcome  is  that  because  of  the 
large  yield,  the  spectrum  is  heavily  shifted  to  lower  frequencies  where  the 
reflection  coefficient  is  rapidly  changing  from  about  1.0  to  some  smaller 
aUssolute  value,  representing  a  failure  of  the  model.  This  matter  is  discussed 
in  more  detail  in  the  following  sections. 


Bakun,  W.  H.  and  L.  R.  Johnson,  1973,  The  deconvolution  of  teleseismic  P 

waves  from  explosions  MILROW  and  CANNIKIN;  Geophys.  J.  R. ,  M,  321-342. 

Springer,  D.  L. ,  1974,  Secondary  sources  of  seismic  waves  from  underground 
nuclear  explosions;  Bull.  Seism.  Soc.  Am.,  64,  581-594. 
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GREELEY  (Figure  6) 

Analysis  of  GREELEY  gives  an  accurate  result  despite  absence  of  a  promi¬ 
nent  spectral  null.  In  this  case  the  maximum- likelihood  technique  is  clearly 
superior  to  standard  methods.  Note,  however,  that  the  theoretical  maximum- 
likelihood  spectrum  is  broader  than  the  observed  spectrum  and  that  the  absolute 
meucimum  of  the  likelihood  function  occurs  for  the  multiple  and  not  for  pP. 

Note  that  in  CANNIKIN,  MILROW,  and  now  GREELEY  there  has  been  a  small  positive 
multiple  around  0.3  seconds.  Experiments  with  artificial  data  show  that  the 
maximum-likelihood  technique  is  unable  to  distinguish  between  a  precursor  and 
a  later  arrival  which  have  the  same  absolute  delay  relative  to  the  largest 
arrival.  (Eliminating  this  ambiguity  is  a  subject  for  future  research.) 
Therefore,  this  multiple  can  be  interpreted  as  a  precursor,  explaining  the 
small  first  motions  observed  at  most  stations  for  these  and  most  other 
explosions. 

LONGSHOT  (Figure  7) 

For  LONGSHOT  agreement  between  theoretically  calculated  spectriom  and 
observed  spectrum  is  good.  A  precursor  is  evident  and  the  pP  delay  is  in 
good  agreement  with  results  from  earlier  workers  and  with  the  first  null  seen 
in  Tables  I  and  II. 

BUKARA  (Figure  8) 

The  event  at  BUKARA  is  different  from  others  in  this  study  because  B  and 
k^  had  to  be  estimated  from  the  data.  The  standard  granite  potential,  together 
with  the  yield  given  by  Marshall  (1972)  of  approximately  30  kt,  proved  unable 
to  reasonadsly  match  the  observed  spectrum  with  any  of  a  complete  range  of  t* 
values.  Marshall  (1972)  stated  that  this  event  was  detonated  in  clay,  a 
medium  very  much  different  from  granite.  One  might  intuitively  expect  little 
overshoot,  a  long  time-constant,  and  low  attenaution  from  this  shot  as  re- 
corded  at  NPNT.  By  trial  and  error  a  good  fit  was  found  to  the  spectral 
envelope  for  B=0,  k  =  7.3,  euid  t*  =  0.16.  These  parcuneters  yielded  a  good 
match  to  the  complete  modulated  spectrum  yielding  a  delay  time  in  good 
agreement  with  Marshall's  (1972)  deconvolution  results  and  with  simple  first 
null  considerations.  This  event,  with  its  large  delay  and  good  high-frequency 
data  is  easy  to  analyze.  The  estimated  reflection  coefficient  is  -0.75,  the 


largest  ever  reported.  This  coefficient  could  result  from  the  (1)  less  rugged 
surface  topography  than  at  other  deep  shot  sites  such  as  NT5;  and  Amchitka, 

(2)  from  the  fact  that  a  small  shot,  being  overhuried,  has  not  crushed  and 
distorted  as  much  of  the  medium  as  larger  shots,  (3)  that  less  of  the  pP 
energy  has  gone  into  spall  since  the  reflected  wave  of  dilatation  may  bo  too 
weak  to  create  substantial  spall  because  the  surface  is  far  from  the  detonation 
point. 

EDAM  (Figure  9) 

In  this  case  good  agreement  is  also  obtained  between  predicted  and  ob¬ 
served  delay  times.  This  shot  is  interesting  because  it  is  relatively  small, 
detonated  at  Yucca  in  a  very  low  velocity  medium,  and  analyzed  only  at  RKON. 
Approximately  the  same  delay  would  he  obtained  by  first  null  considerations. 

MAST  (Figure  10) 

This  event  shows  a  complicated  pattern  of  relatively  weak  arrivals, 
which,  taken  together,  match  the  observed  spectrum  well.  The  strongest  arrival 
of  opposite  polarity  gives  a  delay  of  0.5  seconds,  in  good  agreement  with  the 
predicted  delay.  However,  another  weaker  "reflection"  with  0.95  delay  exists 
that  is  in  fair  agreement  with  the  delay  which  would  be  predicted  from  first 
null  considerations  of  0.8  seconds.  In  this  case  a  clear  first  null  gives  a 
significantly  incorrect  result. 

MI7.ZEN  (Figure  11) 

Here  is  a  simple  case  where  the  predicted  delay  time  is  in  good  agreement 
with  the  maximum-likelihood  estimate  and  with  first  null  analysis. 

PILEDRIVER  (Figure  12) 

For  this  event  good  agreement  exists  between  the  predicted  and  observed 
pP  delay  time  of  2x (0.46/5. 5)  =  0.17  and  0.15  seconds,  respectively.  In  fact, 
the  maximum  likelihood  estimate  is  in  better  agreement  with  the  predicted 
value  than  the  spall-derived  time  of  0.24  seconds  that  Springer  (1974)  pre¬ 
sented.  Although  in  this  case  no  clear  first  null  exists,  a  suggestion  of 


Springer,  D.  L. ,  1974,  Secondary  sources  of  seismic  waves  from  underground 
nuclear  explosions,  Bull.  Seism.  Soc.  Am.,  64.  581-594. 
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one  is  evident  yielding  an  incorrect  delay  of  0.64  seconds.  In  Figure  13 
are  sunmarized  the  foregoing  results  in  a  comparison  of  known  depths  as 
compared  to  depths  calculated  using  known  velocities  and  the  maximum  like¬ 
lihood  estimates  of  t .  We  note  that  the  the  depth  estimate  using  Springer's 
(1974)  T  estimate  has  the  same  size  error  as  the  depth  estimate  from  this 
paper.  The  medium  velocity  for  BUKARA  is  uncertain. 
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IMPLICATIONS  OF  WWSSN  IjONG-RANGF  RECORDINGS  OF  LARGE-EXPLOSION  P  WAVES 


Fiqure  14a  shows  WWSSN  lonq-period  recordinqs  of  CANNIKIN.  Twice  the 
pealc-to-trouqh  time  can  be  estimated  to  be  3.2  seconds,  with  some  care  on  the 
good  signal-to-noise  traces.  In  Figure  14b,  WWSSN  lonq-period  recordings  of 
P  waves  from  a  nearby  earthquake  show  typical  earthquake  periods  of  10  seconds. 
Note  also  that  although  the  CANNIKIN  m.  is  a  full  magnitude  larger  than  that 

D 

of  the  earthquake,  the  lonq-period  amplitude  is  only  2  to  3  times  larger. 
Together,  these  facts  show  that  proportionally  less  long-period  energy  is 
emitted  from  the  CANNIKIN  epicenter  than  from  the  earthquake.  For  CANNIKIN, 

0.8  Hz  is  the  frequency  where  the  amplitude  spectrum  of  scaled  granite-reduced 
displacement  potential  falls  to  one-half  of  the  low-frequency  limit.  A  peak 
in  the  theoretical  spectrum  occurs  at  about  3-second  period.  The  magni¬ 
tude  6.1  earthquake  might  well  have  a  corner  frequency  lower  than  0.3  to 
0.8  Hz.  If  so,  this  alone  might  account  for  the  difference  in  the  observed 
frequency  content  of  the  signals.  However,  the  qranite-reduced  displacement 
potential  is  not  compatible  with  observed  CANNIKIN  waveforms,  unless  the  pP 
phase  cancels  much  of  the  low-frequency  signal.  Without  pP,  the  explosion 
waveform  would  be  predicted  to  show  periods  similar  to  those  of  the  Amchitka 
earthquake.  Thus,  perhaps  pP  is  basically  responsible  for  the  difference  in 
signal  period  between  CANNIKIN  and  the  earthquake. 

Figures  15a  and  15b  show  signals  from  JORUM  at  NTS  and  from  the  nearby 
San  Fernando  earthquake  of  comparable  m^^.  Again,  the  long-period  P-wave 
signals  for  the  explosion  are  of  higher  frequency  and  are  10  to  20  smaller  in 
amplitude.  The  granite  potential  for  JORUM  places  the  one-half  amplitude 
point  for  the  spectrum  at  1.4  Hz,  and  the  spectral  peak  at  1.7  sec.  According 
to  Wyss  and  Hanks  (1972),  the  corner  frequency  for  the  San  Fernando  earthquake 
is  approximately  0.15  Hz,  again  showing  that  the  difference  in  source- spectra 
alone  is  sufficient  to  qualitatively  explain  the  difference  between  the  long- 
period  signals.  However,  quantitatively  we  shall  see  thc.t  we  require  the 
effect  of  pP  to  reduce  the  theoretical  period  expected  from  underground 


Wyss,  M.  emd  T.  C.  Hanks,  1972.  The  source  parameters  of  the  San  Fernando 
earthquake  Inferred  from  teleseismic  body  waves.  Bull.  Seism.  Soc.  Am., 
62,  591-602. 
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Figure  14a.  WWSSN  long-period  recordings  of  P  from  CANNIKIN.  Gains 

as  printed  are  1.63  times  those  indicated  on  the  figure  which 
are  those  given  on  the  original  film. 
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explosions  to  be  in  aqreement  with  observed  values. 

Fiqures  lfia,b,c  show  theoretical  W1\’SSN  lonq-period  waveforms  for  a  ranee 
of  yields  from  20  to  5000  kt  and  for  t*  =  .25,  and  .5.  The  delay  times  for  pP 

1/3 

are  qiven  by  formula  t  *  t seconds  with  Y  in  kilotons.  Calculations  are 
qiven  for  no  pP  and  with  =  0,03,  0.05,  O.l,  0.15,  0.2.  =  0.03  corrt'- 

s{xjnds  to  an  underburied  shot  in  qranite  (a  shallower  PlbKDRIVER) ,  0.20  could 
correspond  to  shots  in  some  of  the  very  low  velocity  sediments  at  Yucca  Plats, 
such  as  EDAM  or  MIZZEN.  The  waveforms  are  calculated  usinq  the  subroutine  PMPHL 
in  programs  WHAMPFS.  The  technique  entails  computing  the  complete  theoretical 
amplitude  spectrum  as  the  product  of  the  von  Seggern-Blandford  source  spectrum, 
an  attenuation  operator  exp  (-t.!t*/2)  and  the  WWSSN  lonq-period  instrument  re- 
sixjnse.  The  phase  is  calculated  as  the  Hilbert  transform  of  the  log  spectrum 
and  the  complete  spectrum  inverted  to  the  time  domain.  The  result  is  convolved 
in  the  time  domain  with  delta  functions  to  give  the  direct  wave  and  pP  as 
appropriate.  For  Figures  16a, b,c  a  reflection  coefficient  of  -1.0  is  assumed. 
Note  that  in  Figure  16  this  technique  automatically  yields  an  increasing  time 
delay  for  first  motion  (as  appropriate  from  causality  considerations'  as  the 
yield  and  t*  becomes  larger. 

To  the  left  of  each  waveform  in  Figure  16  is  plotted  numlM^rs  for  param¬ 
eters  derived  by  automatic  processing.  The  period,  T,  is  taken  to  be  twice 
the  time  interval  between  the  first  peak  and  subsequent  trough.  In  the  case 
of  no  pP  this  assumption  probably  results  in  a  period  which  is  qualitatively 
too  long.  For  the  logarithm  of  the  amplitude,  log  (A),  the  amplitude  is  taken 
from  the  sanK"*  peak-to-trough  excursion.  For  log  (A/T)  ,  A  is  corrected  for 
the  WWSSN  response  at  period  T.  The  integral  from  0.5  to  2.0  Hz  is  calculated 
by  dividing  the  calculated  amplitude  spectrum  by  the  WWSSN  long-period  system 
response,  squaring  the  result,  and  integrating  from  0.5  to  2.0  Hz.  The  result 
is  the  root-mt?an-square  true  ground  displacement  from  0.5  to  2.0  Hz.  The  numbi-i 
plotted  is  the  logarithm  of  this  number. 

Appropriate  waveforms  for  events  in  this  study  are  indicated  by  event  nam<’s 
next  to  the  waveform.  PILEDRIVER  appears  on  two  waveforms  tx'cause  its  scaled 
delay  time  lies  midway  between  .03  and  .05.  Note  also  that  BOXCAR,  HENHAM, 
GREEl.EY,  dORHM,  and  HANDLEY  all  have  similar  source  paranu'ters  and  fall  on 
the  same  waveform. 
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CANNIKIN  provides  the  best  data  available  for  a  comparison  between  theory 
and  experiment.  In  Figure  16b,  the  predicted  period  is  3.2  seconds,  which  is 
in  perfect  agreement  with  the  period  measured  from  the  WWSSN  film.  Note  also 
that  the  detailed  shape  of  the  observed  waveform  agrees  with  the  theoretical 
prediction.  It  does  not  agree  with  the  theoretical  prediction  for  no  pP  in 
Figure  16a. 


Figures  17a, b,c  explore  the  results  of  setting  B=0  in  the  reduced  dis¬ 
placement  potential.  In  Figure  17b  the  period  for  CANNIKIN  becomes  too  long, 
and  the  trough  too  shallow  to  match  observations.  Similar  results  are  obtained 

for  other  values  of  k  above  the  granite  value  of  16.fi.  A  non-zero  value  for 

o 

B,  implying  an  overshoot  in  the  RDP,  is  necessary  to  match  observations.  Fig¬ 
ures  18a, b,c  give  the  results  if  the  pP  reflection  nas  an  eunplitude  of  only 
-0.5  with  respect  to  the  direct  wave.  Again,  for  CANNIKIN  in  Figure  18fa,  the 
period  is  too  long  and  the  trough  has  too  low  an  amplitude.  However,  earlier 
in  this  study,  and  in  other  studies  0.5  has  been  shown  to  be  a  reasonable 
estimate  for  a  reflection  coefficient  at  frequencies  above  0.5  to  1.0  Hz. 


This  paradox  might  be  resolved  by  making  an  analogy  with  results  in 
acoustics  for  plane  wave  reflection  from  rough  surfaces  as  a  function  of  fre¬ 
quency  (Tolstoy  and  Clay  (1966).  Using  an  analytical  expression  of  the  form 

found  in  Tolstoy  and  Clay,  assume  the  reflected  pulse  differs  from  the  direct 

2  2 

pulse  by  the  spectral  amplitude  transfer  function  0.5exp(-f  /f^)  +  0.5.  Note 

this  implies  perfect  reflection  at  f  =  0,  and  a  reflection  coefficient  of 

0.5(e  ^  +  1)  =  0.68  at  f  =  f  .  The  transition  is  rapid  from  1.0  to  0.5  around 

r 

f  =  f^.  Frasier  (1972)  gave  evidence  that  pP  has  proportionally  more  low- 
frequency  amplitude  than  P. 


To  calculate  the  complete  waveform  resulting  from  the  lower  frequency 
pP,  the  previously  used  pP  delta  function  was  replaced  with  the  minimum  phase 
wavelet  with  the  spectrum  shown  above.  Then,  as  before,  the  initial  "P"  delta 


Tolstoy,  I.,  and  C.  I,  Clay,  1966,  Ocean  Acoustics.  Theory  and  Experiment 
in  underwater  sound:  New  York,  McGraw-Hill. 


Frasier,  C.  W. ,  1972,  Observations  of  pP  in  the  short-period  phases  of  NTS 
explosions  recorded  at  Norway;  Geophys.  J.  R.  Astr.  Soc. ,  31,  99-109, 
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Figure  17a.  WWSSN  long-period  theoretical  waveforms  for  a  granite  reduced 

displacement  potential  modified  by  setting  B=0»  assuming  perfect 
reflection  for  pP. 
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Figure  18c.  WWSSN  long-period  theoretical  waveforms  for  a  granite  reduced 
displacement  potential.  pP  reflection  coefficient  is  -0.5. 
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function  plus  this  minimum  phase  wavelet  appropriately  delayed  were  convolved 
with  the  overall  source  mechanism-absorption-instrianent  response  wavelet. 


The  results  are  seen  in  Figures  19a, b,c  and  20a, b,c  for  =  1.0  and 
0,25  respectively.  For  CANNIKIN,  perfect  reflection  existing  only  below 
f  =  1.0  and  perhaps  only  below  0.25  gives  satisfactory  agreement  with  ob¬ 
servations.  Even  with  f^  =  0.25,  agreement  is  much  better  with  Figure  16 
where  the  reflection  coefficient  is  -1.0  than  with  Figure  18  where  it  is  -.5. 
This  analysis  strongly  suggests  that  the  reflection  coefficient  for  pP  from 
CANNIKIN  is  close  to  1.0  for  frequencies  somewhere  below  1.0  Hz. 

Adequately  heuid-digitizing  high-frequency  signals  from  explosions  on 
long-period  records  is  difficult.  For  this  reason,  attempts  were  made  to 
transform  the  LRSM  signals  of  explosions  studied  in  this  report  into  the 
form  they  would  have  had  if  they  were  recorded  on  a  WWSSN  long-period 
system.  To  do  this  we  divided  the  complex  observed  signal  spectrum  by  the 
complex  LRSM  response,  and  then  multiplied  it  by  the  complex  WWSSN  response. 
In  practice,  the  system  phase  responses  are  determined  from  their  log  am¬ 
plitude  responses  by  means  of  the  discrete  Hilbert  transform.  This  process 
is  applied  only  to  the  vertical  components,  and  only  at  those  stations  for 
which,  after  deconvolution,  there  seems  to  be  a  detectable  long-period 
signal.  The  signals  were  aligned  using  the  start  of  the  short-period 
signal.  The  summed  results  are  shown  in  Figure  21.  Periods  of  resulting 
signals  have  been  measured  euid  are  compared  to  the  theoretical  values  in 
Tedile  III,  which  also  includes  a  small  plot  of  the  period  data.  While  the 
agreement  is  good,  there  seems  to  be  a  general  bias:  the  observed  periods 
after  deconvolution  are  shorter  than  those  theoretically  predicted. 

Significantly,  in  this  respect,  the  deconvolved  period  for  CANNIKIN  is 
2,5  seconds,  compared  to  the  3.2  seconds,  observed  directly  from  the  film 
»diich  agreed  with  theory.  Perhaps  there  is  some  effect  of  the  large  long- 
period  noise  level  (most  of  which  is  probably  LRSM  system  noise)  in  the 
deconvolved  data  that  shortens  the  observed  period.  Experiments  using 
the  high  dynamic  range  SRO  stations  may  resolve  this  question. 

Contours  of  log (A)  and  T  as  function  of  yield  and  are  drawn  in 
Figure  22  for  t*  »  0.25  sec  <uid  0.5  sec.  The  values  contoured  are  taken  from 
Figure  16,  the  granite  parameters  with  perfect  reflection.  Note  that  for 
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:ical  waveforms  for  a  granite  reduced 


pP_reflection  coefficient  varies 
25^)  +0.5. 


Figure  21.  LRSM  short-period  vertical  recordings  transformed  to  the 
appearance  they  would  have  on  WWSSN  long-period  systems, 
aligned  on  the  short-period  arrival  times  and  sumned  over 
the  network.  One  sum  trace  is  presented  for  each  event 
in  the  report. 
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fixed  t*  the  families  of  contours  are  generally  not  parallel  so  that  specifi¬ 
cation  of  an  amplitude  and  a  period  will  determine  yield,  and  t  and,  there- 

o 

fore,  the  pP  delay.  The  amplitude  graph  for  t*  *  0.25  could  be  calibrated  by 
CANNIKIN,  and  for  t*  *  0.5,  JORUM  or  some  other  large  NTS  shot.  Clearly, 
however,  the  reduced  displacement  potential  must  be  known  a-priori  to  be  sat¬ 
isfactory  for  the  test  site  under  consideration. 

This  technique  would  be  difficult  to  apply  at  present  because  of  the 
poor  time  resolution  and  signal- to-noise  ratio  on  WWSSN  long-period  traces. 
These  problems  may  possibly  be  overcome  to  some  degree  by  deconvolved  SRO 
data,  and  by  arrays  to  which  multi-channel  filtering  may  be  applied  to  reduce 
the  low-frequency  noise. 


CONCLUSIONS  AND  SUGGESTIONS  FOR  FURTHER  RESEARCH 


In  conclusion,  if  the  parameters  of  the  source  signal  spectrum  are  speci¬ 
fied  correctly,  then  the  time  delays  and  amplitudes  for  a  single  reflection 
can  be  estimated  by  maximum  likelihood.  Difficulties  exist  in  specifying  the 
initial  values  for  the  parameters  B,  C  and  t*.  However,  this  could  be  accom¬ 
plished  by  searching  the  likelihood  over  B,  C  and  t*  but  in  results  not  re¬ 
ported  here  this  has  yielded  unreasonable  values  for  these  parameters.  The 
problem  of  estimating  the  source  spectrum  for  BUKARA  was,  however,  satisfac¬ 
torily  solved  in  this  report  by  an  heuristic  iterative  approach. 

Another  potential  problem  is  the  presence  of  multiple  arrivals  when 
the  model  is  essentially  a  single  reflection  spectrum.  In  the  simulated 
example  presented  here  this  seemed  to  be  a  less  important  problem  as  the 
single  reflection  likelihood  function  successfully  resolved  a  double 
reflection  example. 

Nonetheless  perhaps  another  approach  to  the  multiple  arrival  problem 
should  be  tried.  Conceivably,  this  was  the  difficulty  with  CANNIKIN  because 
five  separate  arrivals  of  high  likelihood  were  detected  and  the  single  arri¬ 
val  model  may  be  inadequate  in  so  complex  a  situation.  Research  also  seems 
warranted  on  the  problem  of  distinguishing  precursors  from  later  multipaths. 

Another  source  of  difficulties  with  CANNIKIN  could  be  that  the  pP 
reflection  coefficient  is  varying  with  frequency  through  the  dominant  part 
of  the  signal  spectrum.  This  variation  violates  the  hypothesized  model  of 
a  constant  reflection  coefficient.  Research  could  be  carried  out  with  a  more 
satisfactory  model  to  determine  if  this  problem  with  CANNIKIN  can  be  solved. 

Various  informal  attempts  have  been  made  to  estimate  the  stability  of 
the  results  to  the  a-prlori  estimates  of  B,  C,  and  t*.  In  general  variations 
of  20  to  30%  in  these  parameters  do  not  significantly  change  the  final 
results.  Still,  however,  a  more  complete  investigation  of  this  question 
seems  warranted. 

Finally,  the  most  promising  result  in  this  study  is  the  successful 
determination  of  the  delay  for  PILEDRIVER.  This  effort  is  of  great  impor¬ 
tance  because  many  of  the  explosions  in  the  USSR  are  apparently  detonated 
in  high-velocity  media,  and  detonated  at  depths  where  the  delay  time  is  very 
short.  The  maximum-likelihood  technique  should,  therefore,  be  tried  out  on 
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